#' This function performs differetial analysis by fitting read
#' Perform likelihood ratio tests and extract the differential
#' analysis results
#'
#' This function performs likelihood ratio tests for given
#' coefficinets contrasts after fitting read counts to GLM by
#' \code{\link{DBanalysis}}. \code{DBresult} extracts the
#' diffential analysis results of given contrasts for all
#' genomic features or genomic features with significant
#' differential events. \code{DBresult.cluster} returns similar
#' results while the results only contain genomic features belong
#' to a given cluster.
#'
#' @name DBresult
#'
#' @param object a \code{TCA} object, for \code{DBresult},
#' \code{DBanalysis} should already be called on the object;
#' for \code{DBresult.cluster}, both \code{DBanalysis} and
#' \code{timeclust} should be already called.
#'
#' @param group1 character string giving the level to be compared,
#' that is the denominator in the fold changes.
#'
#' @param group2 a character vetor giving other levels to compared
#' with \code{group1}. that are numerator in the fold changes.
#'
#' @param contrasts a character vector, each charcter string in
#' the vector gives a contrast of two groups with the format
#' group2vsgroup1', group1 is the denominator level in the fold
#' changes and group2 is the numerator
#' level in the fold changes.
#'
#' @param p.adjust character string specifying a correction method
#' for p-values. Options are 'holm', hochberg', 'hommel',
#' 'bonferroni', BH', 'BY', 'fdr', 'none'.
#'
#' @param top.sig logical if TRUE, only genomic regions with
#' significant differential events will are returned. Significant
#' differential events are defined by log2-fold changes,p-values or
#' adjusted p-values.
#'
#' @param pvalue character string specify the type of p-values
#' used to define significant differential events('\code{PValue}'
#' or adjusted p-value 'paj')
#'
#' @param pvalue.threshold a numeric value giving threshold of
#' selected p-value, Significant differential events have lower
#' (ajusted) p-values than the threshold.
#'
#' @param abs.fold a numeric value, the least absolute log2-fold
#' changes
#'
#' @param direction character string specify the direction of fold
#' changes ('\code{up}' (positive fold changes), \code{down}'
#' (negative fold changes), \code{both}'(both positive and
#' negative fold changes)). Significant events have log2-fold
#' changes exceeding \code{abs.fold} in defined directions.
#'
#' @param cluster an integer, the result tables of genomic features
#' belong to the \code{cluster} are extracted.
#'
#' @param  cmthreshold a numeric value, this argument is applicable
#' only if \code{cmeans}' clustering method is selected when calling
#' \code{\link{timeclust}} function. if not NULL, the result table of
#' genomic features that belong to the defined \code{cluster} and
#' the membership values to this cluster exceed \code{cmthreshold}
#' are extracted.
#'
#' @param result.type character string giving the data type of return
#' value. Options are "GRangesList" and "list".
#'
#' @details This function uses \code{\link{glmLRT}} from edgeR which
#' perform likelihood ratio tests for testing significance of changes.
#' For more deatils,
#' see \code{\link{glmLRT}}
#'
#' @note If not NULL \code{group1}, \code{group2} and \code{contrasts},
#' result tables are extracted from comparisons in \code{constrasts}.
#'
#' @return
#' A list or a GRangesList.
#' If \code{result.type} is "GRangesList", a GRangesList is returned
#' containing the differential analysis results for all provided contrasts.
#' Each GRanges object of the list is one contrast, the analysis results
#' are contained in 4 metadata columns:
#'
#' @return \code{logFC} log2-fold changes of differential event between
#' two tested.
#'
#' @return \code{PValue} p-values.
#'
#' @return \code{paj} adjusted p-values
#'
#' @return \code{id} genomic feature name
#'
#' If \code{result.type} is "list", a List of data frames is returned.
#' Each data frame is one contrast and contains the following columns:
#'
#' @return \code{logFC} log2-fold changes of differential event between
#' two tested.
#'
#' @return \code{PValue} p-values.
#'
#' @return \code{paj} adjusted p-values
#'
#' @return \code{chr}  name of the chromosomes
#'
#' @return \code{start} starting position of the feature in the
#' chromosome
#'
#' @return \code{end} ending postition of the feature in the chromosome
#'
#' @return \code{id} genomic feature name
#'
#' @author
#' Mengjun Wu, Lei Gu
#'
#' @seealso
#'
#' \code{\link{glmLRT}}
#'
#' @examples
#' data(tca_ATAC)
#' tca_ATAC <- DBanalysis(tca_ATAC)
#' ### extract differntial analysis of 24h, 72h to 0h
#' # set the contrasts using the 'group1' and 'group2' paramters
#' res1 <- DBresult(tca_ATAC, group1 = '0h', group2 = c('24h', '72h'))
#' # one can get the same result by setting the contrasts using hte 'contrasts' parameter
#' res2 <- DBresult(tca_ATAC, contrasts = c('24hvs0h', '72hvs0h'))
#' # extract significant diffential events
#' res.sig <- DBresult(tca_ATAC, contrasts = c('24hvs0h', '72hvs0h'),
#'                    top.sig = TRUE)
#'
#' # extract differntial analysis of 24h, 72h to 0h of a given cluster
#' tca_ATAC <- timecourseTable(tca_ATAC, filter = TRUE)
#' tca_ATAC <- timeclust(tca_ATAC, algo = 'cm', k = 6)
#' res_cluster1 <- DBresult.cluster(tca_ATAC, group1 = '0h',
#'                                  group2 = c('24h', '72h'),
#'                                  cluster = 1)
#'
#'
#'
#' @export
DBresult <- function(object, group1 = NULL, group2 = NULL,
                     contrasts = NULL, p.adjust = "fdr",
                     top.sig = FALSE, pvalue = "paj",
                     pvalue.threshold = 0.05, abs.fold = 2,
                     direction = "both", result.type = "GRangesList") {
  if (is.null(group1) && is.null(group2) && is.null(contrasts)) {
    stop("Either information of groups to compare or \"contrasts\" should be provided")
  }
  if (is.null(contrasts)) {
    if (is.null(group1) || is.null(group2)) {
      stop("One group is missing, two groups should be provided for differential event analysis.")
    }
    if (sum(group1 %in% group2) > 0) {
      warning("Members in group1 are also found in group2, overlapped members are removed from group2.")
      group2 <- group2[-which(group2 %in% group1)]
    }
    contrasts <- contrastNames(group1, group2)
  }
  if (!p.adjust %in% c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")) {
    stop("Method for adjusting P-values should be one of following methods: 'holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none'. Character string is case sensitive.")
  }
  fit <- object@DBfit
  contrast.table <- object@contrasts
  gi <- object@genomicFeature[object@genomicFeature$id %in%
                                row.names(fit$coefficients), ]
  gi <- gi[, c("chr", "start", "end", "id")]
  res <- list()
  for (i in contrasts) {
    tmp <- glmLRT(fit, contrast = contrast.table[, i])
    restmp <- tmp$table[, c(1, 4)]
    adjustp <- p.adjust(restmp$PValue, method = p.adjust)
    restmp <- cbind(restmp, adjustp)
    colnames(restmp)[length(restmp[1, ])] <- "paj"
    restmp <- cbind(restmp, gi, stringsAsFactors = FALSE)
    res[[i]] <- restmp
  }
  if (top.sig) {
    res <- DBresult.filter(x = res, pvalue = pvalue,
                           pvalue.threshold = pvalue.threshold,
                           abs.fold = abs.fold,
                           direction = direction)
  }
  if (tolower(result.type) == "grangeslist") {
    gr <- as(do.call(rbind, unname(res)), "GRanges")
    res <- suppressWarnings(split(gr, rep(names(res), lengths(res))))
  }

  res
}

#' @rdname DBresult
#' @export
DBresult.cluster <- function(object, group1 = NULL, group2 = NULL,
                             contrasts = NULL, p.adjust = "fdr",
                             top.sig = FALSE, pvalue = "paj",
                             pvalue.threshold = 0.05, abs.fold = 2,
                             direction = "both",cluster, cmthreshold = NULL,
                             result.type = "GRangesList") {
  if (length(object@clusterRes@cluster) == 0) {
    stop("No cluster information provided, clustering analysis must be performed first")
  }
  DBres <- DBresult(object, group1 = group1, group2 = group2,
                    contrasts = contrasts, p.adjust = p.adjust,
                    top.sig = top.sig, pvalue = pvalue,
                    pvalue.threshold = pvalue.threshold,
                    abs.fold = abs.fold,
                    direction = direction, result.type = "list")
  names <- names(object@clusterRes@cluster)
  res <- list()
  contrast_name <- names(DBres)
  counter <- 0
  for (i in DBres) {
    restmp <- i
    counter <- counter + 1
    clusters <- object@clusterRes@cluster
    clusternames <- names[which(clusters == cluster)]
    if (!is.null(cmthreshold)) {
      membership <- object@clusterRes@membership[clusters == cluster,
                                                 cluster]
      if (is.null(membership)) {
        stop("No membership matrix found. To get membership matrix, please choose 'cmeans' clustering method when calling timeclust")
      } else {
        clusternames <- clusternames[which(membership > cmthreshold)]
      }
    }
    restmp <- restmp[clusternames, ]
    contrast <- contrast_name[counter]
    res[[contrast]] <- restmp
  }
  if (tolower(result.type) == "grangeslist") {
    gr <- as(do.call(rbind, unname(res)), "GRanges")
    res <- suppressWarnings(split(gr, rep(names(res), lengths(res))))
  }
  res
}

# contrast contrast by given strings, group1 is a string, group2 can be a string or a vector of strings
contrastNames <- function(group1, group2) {
  b <- length(group2)
  name <- vector(mode = "character", length = b)
  for (i in seq_len(b)) {
    n <- paste0(group2[i], "vs", group1)
    name[i] <- n
  }
  name
}

DBresult.filter <- function(x, pvalue = "paj", pvalue.threshold = 0.05,
                            abs.fold = 2, direction = "both") {
  if (abs.fold < 0) {
    err <- paste0("\"abs.fold\" should be postive number.")
    stop(err)
  }
  d <- x
  for (i in seq_len(length(d))) {
    dt <- d[[i]]
    if (direction == "up") {
      dt <- dt[which(dt$logFC >= abs.fold), ]
    }

    if (direction == "down") {
      dt <- dt[which(dt$logFC <= -abs.fold), ]
    }
    if (direction == "both") {
      if (abs.fold == 0) {
        dt <- rbind(dt[which(dt$logFC >= abs.fold), ],
                    dt[which(dt$logFC < -abs.fold), ])
      } else {
        dt <- rbind(dt[which(dt$logFC >= abs.fold), ],
                    dt[which(dt$logFC <= -abs.fold), ])
      }
    }
    dt <- dt[which(dt[, pvalue] < pvalue.threshold), ]
    d[[i]] <- dt
  }
d
}
